### Alizade, Dancygier, Ditlmann
### "National Penalties Reversed"
### Replication Code 
### Figure A3
### For questions, contact jalizade@princeton.edu

# setup
rm(list = ls())
setwd("C:/Users/Jey/Dropbox/WZB/NaturalizationExperiment/Submission/JOP/replication_JOP/data")
library(foreign)
library(lmtest)
library(sandwich)
library(ggplot2)
dat <- read.dta("data_experimental.dta")

# define outcomes
outcomes_1 <- c("e1_legal_dummy", "e1_help_dummy", "e1_posaff_dummy", "e1_legal_sum", "e1_help_sum", "e1_posaff_sum")
outcomes_2 <- c("e2_legal_dummy", "e2_help_dummy", "e2_posaff_dummy", "e2_legal_sum", "e2_help_sum", "e2_posaff_sum")

# convert treatment variable from second experiment to factor
dat$e2_treat[dat$e2_treat=="NA"] <- NA
dat$e2_treat <- factor(dat$e2_treat, levels=c("Vote Intention (Canadian)", "Control (Turkish)", "Vote Intention (Turkish)", "Integration Problems (Turkish)"))


### Figure A3a ###

# regressions
mods_study1 <- lapply(outcomes_1, function(x) {
    mod <- lm(dat[,x] ~ e1_treat_turkish + e1_treat_dual + e1_treat_benefits, data=dat)
    mod_rob <- coeftest(mod, vcov = vcovHC(mod, "HC1"))         
})

# create data frame with estimates
ests_study1 <- data.frame(do.call(rbind, mods_study1)[-seq(1, 21, 4),1:2])
names(ests_study1) <- c("coef", "se")
ests_study1$treat <- rep(c("Turkish", "Dual", "Benefits"), 6)
ests_study1$outcome <- rep(c("Information (binary)", "Help (binary)", "Positive Affect (binary)", "Information (count)", "Help (count)", "Positive Affect (count)"), each=3)
ests_study1$outcome <- factor(ests_study1$outcome, levels=c("Information (binary)", "Help (binary)", "Positive Affect (binary)", "Information (count)", "Help (count)", "Positive Affect (count)"))

# plot
fig_study1 <- ggplot(data=ests_study1, aes(x=treat, y=coef))
fig_study1 <- fig_study1 + geom_point(size = 2.5)
fig_study1 <- fig_study1 + geom_errorbar(aes(ymin=coef-se, ymax=coef+se), size = 1.2, width = 0)
fig_study1 <- fig_study1 + geom_errorbar(aes(ymin=coef-se*1.96, ymax=coef+se*1.96), size = 0.4, width = 0.2)
fig_study1 <- fig_study1 + xlab("") + ylab("Treatment Effect")
fig_study1 <- fig_study1 + geom_hline(yintercept = 0, linetype = "dashed")
fig_study1 <- fig_study1 + scale_y_continuous(breaks = seq(-.1, .2, .1), labels = c("-0.1", "0.0", "0.1", "0.2"), limits = c(-.15, .2))
fig_study1 <- fig_study1 + facet_wrap(~outcome)
fig_study1 <- fig_study1 + coord_flip()
fig_study1


### Figure A3b ###

# relevel treatment variable
dat$e2_treat <- relevel(dat$e2_treat, ref="Vote Intention (Turkish)")

# regressions
mods_study2 <- lapply(outcomes_2, function(x) {
  mod <- lm(dat[,x] ~ e2_treat, data=dat)
  mod_rob <- coeftest(mod, vcov = vcovHC(mod, "HC1"))         
})

# create data frame with estimates
ests_study2 <- data.frame(do.call(rbind, mods_study2)[-seq(1, 21, 4),1:2])
names(ests_study2) <- c("coef", "se")
ests_study2$treat <- rep(c("Vote Intention\n(Canadian)", "Control\n(Turkish)", "Integration\nProblems (Turkish)"), 6)
ests_study2$treat <- factor(ests_study2$treat, levels=c("Integration\nProblems (Turkish)", "Vote Intention\n(Canadian)", "Control\n(Turkish)"))
ests_study2$outcome <- rep(c("Information (binary)", "Help (binary)", "Positive Affect (binary)", "Information (count)", "Help (count)", "Positive Affect (count)"), each=3)
ests_study2$outcome <- factor(ests_study2$outcome, levels=c("Information (binary)", "Help (binary)", "Positive Affect (binary)", "Information (count)", "Help (count)", "Positive Affect (count)"))

# plot
fig_study2 <- ggplot(data=ests_study2, aes(x=treat, y=coef))
fig_study2 <- fig_study2 + geom_point(size = 2.5)
fig_study2 <- fig_study2 + geom_errorbar(aes(ymin=coef-se, ymax=coef+se), size = 1.2, width = 0)
fig_study2 <- fig_study2 + geom_errorbar(aes(ymin=coef-se*1.96, ymax=coef+se*1.96), size = 0.4, width = 0.2)
fig_study2 <- fig_study2 + xlab("") + ylab("Treatment Effect")
fig_study2 <- fig_study2 + geom_hline(yintercept = 0, linetype = "dashed")
fig_study2 <- fig_study2 + scale_y_continuous(breaks = seq(-.1, .2, .1), labels = c("-0.1", "0.0", "0.1", "0.2"), limits = c(-.15, .2))
fig_study2 <- fig_study2 + facet_wrap(~outcome)
fig_study2 <- fig_study2 + coord_flip()
fig_study2
